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Summary 


During the first part of this fourth reporting period research has concentrated on performing a 
detailed evaluation, to zero order, of the guidance algorithm developed in the first period taking the 
numerical approach developed in the third period. A zero order matched asymptotic expansion 
(MAE) solution that closely satisfies a set of 6 implicit equations in 6 unknowns to an accuracy of 
10-10, was evaluated. Guidance law implementation entails treating the current state as a new 
initial state and repetitively solving the MAE problem to obtain the feedback controls. A zero order 
guided solution was evaluated and compared with the optimal solution that was obtained by 
numerical methods. Numerical experience shows that the zero order guided solution is close to 
optimal solution, and that the zero order MAE outer solution plays a critical role in accounting for 
the variations in Loh's term near the exit phase of the maneuver. However, the deficiency that 
remains in several of the critical variables indicates the need for a first order correction. During the 
second part of this period, we explored methods for computing a first order correction. 


1. Introduction 


The primary objective of this research effort is the development of near optimal guidance 
solutions for aeroassisted orbit transfer vehicles that are implementable in real time and on-board 
the vehicle. Therefore, solutions that are both near optimal and that require a minimum of 
computation are of primary interest. During the first reporting period, the application of matched 
asymptotic expansions (MAE) to the problem of inclination change with minimum energy loss has 
been explored. A complete set of integrals for the state and costate equations for this problem were 
found 1 ’ 2 . Enforcing the matching conditions, boundary conditions and optimality conditions 
results in a set of 20 nonlinear algebraic equations. The solution of these equations provides the 
information needed to form a guidance algorithm. An attempt to solve this set of equations was 
made during the second reporting period, but was unsuccessful due to its complexity 3 . An 
alternative approach that greatly simplifies the solution procedure was taken in the third period 
during which it was found that, by exploiting the structure of the MAE solution procedure, the 
original problem could be further reduced to a set of 6 implicit equations in 6 unknowns. During 
the first part of the fourth reporting period research has concentrated on performing a detailed 
evaluation, to zero order, of the guidance algorithm developed in the first period taking the 
numerical approach developed in the third period. Guidance law implementation entails treating the 
current state as a new initial state and repetitively solving the MAE problem to obtain the feedback 
controls. The optimal solution, obtained using numerical methods, serves as a reference solution 
to compare with the guided solution. During the second part of the fourth reporting period, we 
have explored methods for obtaining first order correction. 

A straightforward approach entails numerical solution of the first order Euler-Lagrange 
equations (state and costate system). This approach requires that the state transition matrix be 
obtained by quadrature, and a second quadrature is needed to deal with the forcing term. The state 
transition matrix is obtained by computing unit responses and collecting them in a matrix. This 
approach appears to be computationally intensive if it has to be performed on line (although the 
quadratures can be performed in parallel), but it offers the potential of precomputing the zero order 
solution and the quadratures and storing them for later use in flight. An alternative approach 
consists of performing a MAE expansion the Hamilton- Jacobi-Bellman partial differential equation, 
and performing a single quadrature to obtain the first order costates. Although this requires only a 
single quadrature, it must be performed on-line at each control update. It has been shown that 
using this approach, the first order correction for the outer costates do not contribute to the 
composite costates, thus it is sufficient to obtain the first order correction for the inner costates 
alone. The current effort entails using this approach to obtain a numerical evaluation of a first 
order guided solution. 
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2. Current Research Results 


2.1 Zero Order MAE Approximate Solution 

The numerical method used to obtain the zero order MAE solution is developed in Ref. 4. A 
summary of this procedure and the results from Ref. 1 are presented below. A complete set of 
integrals for the state and costate equations for the problem of inclination change with minimum 
energy loss, were found in Ref. 1. Enforcing the matching conditions, boundary conditions and 
optimality conditions results in a set of 20 nonlinear algebraic equations. The solution of these 
equations provides the information needed to form a guidance algorithm. By exploiting the 
structure of the MAE solution procedure the original problem is simplified by further reducing it to 
a set of 6 implicit equations in 6 unknowns. The unknowns are the common parts of the inner and 
outer solution (which are equal to one another). The iteration procedure involves repeated solution 
of the inner and outer problems using the common parts as artificial boundary conditions. The 
common parts are adjusted in the iteration process until the actual boundary conditions are satisfied 
by the composite solution. The matching conditions are enforced at each iteration by simply 
equating the inner and outer common parts. 

Summary of Results from Ref. 1 

The reduced three dimensional point mass equations of motion for a lifting vehicle over a 
spherical non-rotating planet are given by: 

du/dh = -Bu(l+A. 2 )e- h / £ /eE*siny - 2/(l+h) 2 (1) 

dy/dh = BA,cospe _h/e /esiny + [l/(l+h)-l/u(l+h) 2 ]coty (2) 

dvj//dh = BA.sin|ie h/e /esinycosy (3) 

where the parameter e is the ratio of the atmospheric scale height to the minimum trajectory radius 
r s . For Earth the value of £ is approximately 1/900. The definitions for the remaining variables are 
given below. 

The controls are the normalized lift coefficient X and the bank angle (I which are assumed not 
to be beyond there limits. The optimal control is obtained as a function of the state and costate 
variables by solving the optimality conditions. The resulting expressions are: 
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X = E*(p Y cos|i+pysinp/cosy)/2up u 


(4) 


tanp. = py/pycosy (5) 

where p u , Py and Py are the associated costate variables. 

The zero order equations for the outer problem can be obtained by simply taking the limit as 
e approaches zero on the right hand side of Eq's. 1-3. The general solution for the outer system to 
zero order in e is given by: 


u° = 2[ci + l/(l+h)] 

(6) 

cosy 0 = c 2 /(l+h)u ol/2 

(7) 

O 

II 

o 

(8) 

PS = -a2/2u°+ai 

(9) 

P® = a2tany° 

(10) 

Py = 

(ID 

The dimensionless variables h and u appearing above are defined in terms of the original 
dimensional state variables as: 

h = (r-r s )/r s 

(12) 

u = V 2 /g s r s 

(13) 

where r s is the reference radius corresponding to the lowest point in the trajectory. In the inner 
region, where the aerodynamic force is dominant, a new stretched altitude h and state variables w 
and v are defined as: 

w = Be h , h = h/e, B = CLP s s/2m|5 

(14) 

v = E*ln(l/g s r s u), E* = Cl*/Cd* 

(15) 
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The zero order equations for the inner problem can be obtained by expressing the system in 
Eq's. 1-3 in terms of the inner variables appearing in Eq's. 14-15, and simply taking the limit as e 
approaches zero on its right hand side. The general solution for the inner system to zero order in e 
is given by: 


7° = -ki\|/° 2 /2 + k2V°+ k3 

(16) 

= [kjya/6 - k 2 \j/° 2 /2 - k 3 \p+ k 4 ]/o 

(17) 

v° = (a+l/o)\j/°+ o[(y°ki-k2) 3 ]/3ki + ks 

(18) 

P° = (Pw/a)V° + c 

(19) 

Pw = constant 

(20) 

Pv = constant 

(21) 

where 


ki = P&/2a 2 P? , k 2 = -c/2opf 

(22) 

The superscript "o" denotes the fact that these are the zero order solutions in the expansion 
variable, e = l/(3r s . The overbars in the above equations are used to distinguish an inner solution 

variable from an outer solution variable. 

In the above expressions V is the velocity, r is the radius from the earth center, z is the 
altitude above sea level, (3 is the inverse of the scale height associated for an exponential air density 
model, Ps is the reference air density at r = r s , s is the reference area, m is the mass and Cl* and 
Cd* are the lift and drag coefficients corresponding to the maximum lift to drag ratio. The 
subscript s is used to denote the reference radius value, corresponding to the lowest altitude of the 
trajectory. Note that h=y=0 at this reference radius, and that h was used as the independent 
variable for the original problem and for the outer solution, while \p was used as the independent 


variable for the inner solution. 

The inner solution variables w and v are related to the outer solution variables h and u by 
the transformations in Eq's. 14 and 15. The outer solution costate variables can be expressed in 
terms of the corresponding inner solution costate variables using: 
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po = -E*P$/u, 


p° = -2op? 


(23) 


where a is the normalized horizontal lift component given by 

a = [1/(1 + k 2 2 + 2kik 3 )] 1/2 (24) 

A uniformly valid zero order approximation for the optimal solution is constructed by 
adding the outer and inner solutions, and subtracting their common parts 

x c (h) = x°(h) + x°(h) - x°(0) (25) 

The superscript c is used to denote this composite solution, and x is used to denote any of the 
variables of interest. Note that the composite solution must have h as the independent variable, to 
be consistent with the original MAE problem formulation. Also note that the matching conditions 
require x°(0) = x°(°°), and that the boundary conditions are satisfied using the composite solution. 

The entry and exit Keplerian arcs (outer solution) are different in their orbital parameters 
thus there are two separate outer solutions (called the left and right solutions in Ref. 1). In total, 
there are 18 integration constants (twelve from the two outer solutions and six from the inner 
solution) that need to be evaluated to obtain a composite solution. The reference radius (r s ) and the 
horizontal component of lift (o), which was shown to be constant to zero order in e, bring the total 
number of unknowns to 20. The left and right matching conditions and the boundary conditions 
for the composite solution provide 18 algebraic equations. Two additional equations are obtained 
by using H<j = 0 (Eq. 24 above) and the fact that r = r s when 'f = 0. This gives a total of 20 

coupled algebraic relations between the constants and an iterative method is needed to calculate 
them. 

Figure 1 presents a summary of results from Ref 4. The converged zero order MAE 
solution for y, v, Py, Pu and the controls X and p indicate that the composite solutions for y and p Y 
are significantly different from the inner solution. In particular, the major variation in p Y is due to 
the outer solution, which (as explained in Ref. 1) in effect amounts to a correction for the large 
variation in Loh's term during the exit phase. The implications of this important correction have 
been fully addressed in Ref. 1, and will not be elaborated here. The normalized lift control X is 
always near 1 corresponding to flight a near maximum lift to drag ratio. The bank angle |i is 
always near 90 degrees indicating that most of the aerodynamic force is utilized in performing the 
turn. 
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2.2 Zero Order Guided Solution 


Closed loop guided solutions are obtained by using the optimal control expressions given in 
Eq.s 4 and 5. These expressions involve both the states and the costates thus knowledge of both 
states and costates is needed to evaluate the controls. Assuming that the states are available for 
feedback, only estimates of the costates are required at each control computation along the 
trajectory. Feedback implementation entails treating the current state (from the simulation) at each 
control update as a new initial state, and calculating the costate values corresponding to the same 
time instant. The estimate for these costates (to zero order) are obtained by repetitively solving the 
zero order MAE problem. 

At the first step, an initial guess and boundary conditions are supplied to initiate the procedure 
of obtaining a zero order MAE converged solution. Next, the costate expressions in Eq's. 9-11 
and 19-21 can be evaluated as a function of the corresponding independent variables and used in 
Eq. 25 to construct the composite costate expressions which are in turn used in Eq's 4 and 5 to 
compute the controls. When a predetermined time increment has been reached, the current states 
are used as initial conditions for the next MAE calculation. It follows that the initial guess is 
available in every step of the zero order MAE calculation after the first. 

For between-updates control computation, the integration constants from the last update are 
used. The transformations defined in Eq's. 12-15 are used to transform the simulated dimensional 
variables into the inner and outer dependent and independent variables. The transformed altitude h 
is used in Eq's. 6,7,9 and 10 to compute the left and right outer costates. The heading \ji° is used 
in Eq. 19 to evaluate Py, and Eq. 25 is used to calculate the composite costates. The composite 
costates and the current y and u (from Eq. 13) are used in Eq's. 4 and 5 to evaluate the optimal 
controls between the MAE update calculations. 

During the exit phase, the left outer solution is discarded, and matching is required only 
between the right outer solution and the inner solution. In this case, the constants of integration for 
the inner solution are viewed as free parameters used to satisfy the boundary conditions. 

Numerical Results 

Fig's 2 through 6 present a comparison between the optimal solution and the zero order 
guided solution for A\y=20°. The corresponding control time histories are given in Fig's. 7 and 8. 
Loh's term (corresponding to the optimal solution) is given in Fig. 9. The time increment between 
guided solution updates is 5 seconds, with the control computed between updates at every 
integration step following the procedure described earlier. These results indicate that the guided 
solutions and the optimal solutions are in a very good agreement throughout the trajectory. The 
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error that does exist in some of the guided solution variables, for example the y and the altitude (z) 
solutions, indicate the need for a first order correction. Fig’s 4 and 7 clearly indicate how the 
variation in Loh's term near the end of the trajectory is partially accounted for in the zero order 
guided solution. Specifically, the normalized lift coefficient X given in Fig. 7 does not saturate in 
the exit phase as in previous studies 5 , but reduces to near zero. Fig. 10 compares the velocity 
histories near the end of the trajectory. The optimal and guided solution values for terminal 
velocities are 6751 m/sec and 6736 m/sec respectively. The change of the guided value from the 
optimal one is 15 m/sec. These results are to be presented at the AIAA Guidance, Navigation and 
Control Conference, Hilton Head, SC, August 10-12, 1992 6 . 
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2.3 First Order Correction 


The results of Section 2.2 indicate the need for a first order correction due to the deficiency 
that remains in some of the critical variables, for example the y and the altitude (z) solutions (Fig's. 
2 and 6). A straightforward approach to obtain a first order correction is by retaining the 0(e) 
terms in the expansion of equation 1-3 and the corresponding costate equations. This result in a 
nonlinear zero order system the solution to which was obtained analytically in Ref. 1 , and a linear 
system with variable coefficients for the first order terms of the expansion. The solution to the first 
order linear system can be obtained by numerically constructing the state transition matrix and 
evaluating the first order costates at the initial time. The costates values to first order are used in 
the control expressions 4-5 which provide the feedback control functions for the guided solution. 
The drawback of this approach is that the state transition matrix must be computed from unit 
responses. This means that 2n quadratures are required (where n is the state dimension), followed 
by an additional quadrature to handle the forcing term. The 2n quadrature for the state transition 
matrix computation can be performed in parallel. The advantage is that it offers the potential of 
precomputing the zero order solution and the quadratures and storing them for later use in flight 7 . 
This issue will be explored during the next reporting period 

An alternative approach to obtain guided solutions to higher orders is to assume that the states 
are available for feedback, thus only estimates of the costates are needed at each control 
computation along the trajectory. In Ref. 5 a regular perturbation method was used to obtain a first 
order correction for the costates by performing quadratures at each control update. However, the 
problem addressed is actually the inner part of a two time scale analysis based on singular 
perturbation theory (Ref. 1). Thus there is no correction for the variation in Loh's term in the zero 
order solution. In theory, the solution in Ref. 5 is not a uniformly valid approximation to any 
order in e. 

Summary of Results from Ref. 5 

The dynamic system is given by 

y = f(y,u) + eg(y) (26) 

where y is an n-dimensional state vector, u is an m-dimensional control vector, £ is a small 
parameter, x is the independent variable, y = dy/dx, x is the initial state and t is the initial value of 
the independent variable. The optimization problem is to find u to minimize J=<|>(yf) subject to 
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dynamics of Eq' (1) and the terminal constraint \|/(yf)=0. The Hamilton-Jacobi-Bellman partial 


differential equation (HJB-PDE) is 


P t = - H°P l 

(27. a) 

H°P l = min H = P x (f°P t + eg) 

ueU ^ 

(27. b) 


where U is a class of continuous-bounded controls, f°P l = f[y,u°Pt(x,t)], and u 0 P l (x,t) is given by 
the optimality condition H u = 0, assuming that H uu is positive definite. P(x,t) is the optimal return 
function defined as the performance index for an optimal path starting at x and t and satisfying the 
terminal conditions. P(x,t) is expressed as a series expansion in £ as 

oo 

P(x,t) = X Pj (M)e j (28) 

j=0 

and the optimal control as 

oo 

u°P l (x,t) = u°P l (x,P x ,t) = X u j( x >0 £ ^ (29) 

j=o 

where the expression for u°P l (x,t) is obtained by substituting Eq. (28) into the expression obtained 
from the optimality conditions u°P t (x,P x ,t) and expanding the function. Hence, by determining the 
partials Pj x , an optimal control law in a feedback form can be constructed. By expanding f(y,u) 
and substituting the expansions for f(x,t) and P(x,t) into Eq. (27.b), the expansion of the HJB- 
PDE is given by 

oo oo oo 

Pt = X Pj. £j = - ( X Pj. £j )( eg + X f je j ) (3°) 

j=o j=o j=0 

Assuming that an analytic solution for the zero order problem is available, the first order HJB-PDE 
is given by 


Pi, + Pi,fS pt =-Po,g( x ) ( 31 > 

Partial differential equations of this type are solved by the method of characteristics. The 
characteristic curves of the HJB-PDE, for any order of pj, are given by the zero order optimal 
trajectory 
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(32) 


y = f 0 


whose solution is denoted as y 0 (x; x, t). Then the solution for P, in Eq. (31) is given by 


*i(x, 0 =-j Po,g(y)dx 


(33) 


The partials P lx (the costates), which are needed to construct the optimal control u b are given by 
differentiating Eq. (33) with respect to the arbitrary current conditions x as follows: 


Pu=9p 1 /3x = - 


-fa,. 


g (y)]/axdx - [p 0l g(y)] tf at f /ax 


(34) 


and are valid only at the initial time t. 

First Order Correction in the Context of the MAE Problem 

In Ref. 1 it was shown that the problem in [5] is actually the inner part of a two time scale 
analysis based on singular perturbation theory where the constant Loh's term approximation, for 
which complete analytic results are available, corresponds to the zero order solution in [5]. It 
follows that to generalize the approach in [5] to two time scale problems, the analysis in [5] has to 
be carried out for the inner and outer sub-problems separately, and the boundary conditions have to 
be satisfied by the composite solution of the return function. Following the assumption that led to 
Eq. 28 in the context of the MAE problem, each term in the expansion is actually the composite 
solution of the return function to any order in e. The first order term of the composite return 
function is constructed by adding the outer and inner first order solutions and subtracting the 
common part 


Pi(x,t) = P°(x,t) + Pi(x,t/e) - P°(x,0) 


(35) 


In Appendix A it is shown that the first order return function is constant with respect to the 
independent variable, thus it does not contribute to the first order composite return function. 
Hence the inner first order return function, evaluated along the zero order inner analytic solution 
(Ref. 1) and satisfying the boundary conditions, serves as the composite first order correction for 
the return function. Consequently, the analysis in [5] can be used to obtain a first order correction 
for the costates with some modifications that pertain to two time scale type problems. Specifically, 
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Eq. (28) is understood to be the expansion of the inner return function, and the rest of the analysis 
in [5] is valid up to Eq. (33) where the integration is performed along the zero order inner 
characteristic curve. The main conceptual difference is in the way the costates are obtained in Eq. 
(34). In [5], the derivatives in Eq. (34) are of the zero order solution with respect to the true initial 
conditions because the notion of inner and outer solution do not exist in the regular perturbation 
formulation that was used there. In the context of the MAE problem, the derivatives in (34) are 
understood to be of the zero order inner solution with respect to the true initial conditions. In Ref. 
1 numerical evaluation of the zero order solution amounts to a solution of a set of 20 equations for 
18 integration constants and two parameters. This may suggest a practical way to calculate the 
derivatives in (34) of the inner solution with respect to the true initial conditions. 

3(inner solution) _ 3(inner solution) 9(20 parameters) 

3(true initial conditions) 3(20 parameters) 3(true initial conditions) ' 

The details of this method are given in Appendix A. The composite costates to first order, together 
with the measured states, are used in Eq's. 4 and 5 to obtain the control expressions which are 
used in the simulation program. 

2. Future Research 

During the next reporting period we plan to complete the development and evaluation of the 
guidance algorithm including first order corrections. Both approaches outlined in the preceding 
section will be evaluated, including the possibility of precomputing and storing the zero order 
solution and all the required quadratures for the first order solution for the first method discussed. 
This will form an complete evaluation of MAE methods for skip trajectories from the perspective of 
guidance law development. 
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Appendix A 

Auxiliary Derivations 

A.I The first order outer costate 

Equations (1-3) can be written in the following form 

8y/9h = f(y,h) + g(y,u,h/e)/e (A.I) 

where y is an n-dimensional state vector, u is an m-dimensional control vector, e is a small 
parameter and h is the independent variable. The optimization problem is to find u to minimize 
J=<Kyf) subject to dynamics of (A.I) and the terminal constraint t|/(yf)=0. The outer Hamilton- 
Jacobi-Bellman partial differential equation (HJB-PDE) is 


P t = - H«P‘ 

(A.2) 

H°P ( = min H = P x (f + g°P‘/e) 

(A. 3) 


where U is a class of continuous-bounded controls, g°P l = gty.uop^yi.hi)], and u°P l (yi,hi) is given 
by the optimality condition H u = 0, assuming that H uu is positive definite. P(yj,hj) is the outer 
optimal return function defined as the performance index for an optimal path starting at yj and hj 
and satisfying the terminal conditions. P(yi,hi) is expressed as a series expansion in e 

oo 

P(yi,h i ) = XPj(y..h 1 )e j (A- 4 ) 

j=0 

and the outer optimal control as 

oo 

u°Pt(yi,hi) = u°P t (yi,P y ,h) = £ u J (y 1 ,h i ) e j (A.5) 

j=o 

where the expression for u°P ( (y,h) is obtained by substituting (A.4) into the expression obtained 
from the optimality conditions u°P l (y,Py,h) and expanding the function. Hence, by determining 
the partials Pj y , an optimal control law in a feedback form can be constructed. By expanding 
g(y,u,h/e)/e and substituting the expansions for g(y,u,h/e)/e and P(y,h) into Eq. (A. 3), the 
expansion of the outer HJB-PDE is given by 


- 18 - 



p h = X Pj,e j = - ( x p jy £i )( f + X (gj/e)e J ) 

j=0 j=0 j=0 


(A.6) 


By observation of Eq's (1-3) the functions g(y,u,h/e)/e can be decomposed in the following 
manner 


g(y,u,h/e)/e = g 1 (y,u)e- h / £ /e 


(A. 7) 


where gi(y,u) are the parts of g(y,u,h/e)/e that do not contain h or e. It is clear that the exponential 
part in Eq. (A.7) and its derivative, evaluated at e=0, are identically zero, thus the outer HJB-PDE 
are homogeneous to any order in e and the first order outer equation is given by 


Pi h +Pi y f =0 


(A. 8) 


It follows that the first order return function is constant and does not contribute to the first order 
composite return function. 

A.2 Sensitivity of the inner solution to perturbations in the true initial 
conditions 


In [5] the dynamic equations where transformed into the form of Eq. (26) through the use of 
the following nondimansional variables 


v = ln(V2/gr) , w = C*pS/2mP , dt/dz = (VwV 


(A.9) 


where V is the velocity, g is the gravitational acceleration, r is the radius from the earth center, p is 
the density, S is the reference area, m is the mass, (3 is the inverse of the scale height and z replaces 
time (x) as the independent variable. The zero order dynamic equations in [5] are obtained by 
simply setting e=0 in (26) and are identical to the zero order inner equations in the MAE problem. 
The first order inner costates are calculated using Eq. (34) but in the context of the MAE problem. 

(A. 10) 


Pi y =3P 1 /3y = -J^ 3[Po y g(y)]/9yi dz - [P^gty^zf/dyi 


where the integrand in (A. 10) is the partial derivative of the function Po y g(y) of the zero order inner 
solution, with respect to the true initial conditions yi- Since the true initial conditions are not the 
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same as the inner initial conditions, it is necessary to calculate the derivative of the function Po y g(y) 
in two steps. We rely on the fact that the zero order MAE solution is defined by a set of 20 
parameters that are evaluated so as the true boundary conditions are satisfied by the composite 
solution. This suggests that we first take the derivative of Po,g(y) with respect to this set of 20 

parameters, and multiply it by the sensitivity of these parameters to perturbations in the true initial 
conditions. 

The explicit form of the function Po y g(y). denoted by R t , is given in Ref. [5] as 


R, =|^[Pv(l-2e- 7 )Y+Pg(l-e- 7 )] (A.ll) 

The zero order inner solution for the MAE problem was obtained in Ref. 1 and is summarized in 
Eq's. (16-22). Out of the 20 parameters that define the composite zero order MAE solution in Ref. 
1, only 8 appear in the inner solution (16-22) and in R t . These are k 3 , k 4 ,k 5 , c, P w , P v , O and r s . 
What remains is to take the partial derivative of Rj with respect to each of these parameters and 
multiply the result by the sensitivity of the parameters to perturbations in the true initial conditions. 
Since the velocity V was normalized differently in [5] and in [1], it is necessary to relate between v 
defined in Eq. (15) and v as defined in Eq. (A.9), and between the corresponding costates P v and 
P v . These relations are given by 

v = -v/E* -ln(gr) , P v = -E* P v ( A - 1 2 ) 

Using Eq’s. (16-22) in (A.ll, A.12) gives the explicit dependency of Ri on the above 8 
parameters. A relation between the differentials of w and r as obtained by the use of the definitions 
in (12) and (14), is given by 

dr = -dw/{3w (A. 13) 

and the partials of R, with respect to each of the 8 parameters can now be taken explicitly. 

Finally, it is left to determine the sensitivity of the set of parameters to perturbations in the true 
initial conditions. Denoting the set of parameters by a vector x, and the true initial conditions by a 
vector yj, the set of 20 equations for x can be written in the following form (Ref. 1) 

F(x,y0 = 0 (A. 14) 

Differentiating (A. 14) we have 
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dx/dyi = -tfF/dxl-^F/ayi 


(A. 15) 


where the Jacobian matrix dF/dx is numerically calculated once every control update using the 
converged zero order MAE solution. Finally, the partials of R t with respect to the true initial 
conditions are given by 

3R 1 /9yi = 3R[/8x x dx/dyi (A. 1 6) 

This expression is used in (A. 10) to evaluate the first order correction of the inner costates which 
are used in (4,5) to find the correction for the controls. 
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